PROGRAM ALGO836_TEST
    IMPLICIT NONE
    ! These should be included as a subroutin
    INTEGER, PARAMETER :: COLAMD_KNOBS = 20
    INTEGER, PARAMETER :: COLAMD_STATS = 20

    ! The following is problem specific
    INTEGER, PARAMETER :: A_NNZ  = 11
    INTEGER, PARAMETER :: A_NROW = 5
    INTEGER, PARAMETER :: A_NCOL = 4
    INTEGER, PARAMETER :: B_NNZ  = 4
    INTEGER, PARAMETER :: B_N    = 5
    INTEGER :: ALEN
    INTEGER :: row, col, pp, length, ok
    INTEGER, ALLOCATABLE, DIMENSION(:) :: A, B, p,q, perm, stats
    DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: knobs
    !EXTERNAL colamd_recommended
    INTERFACE
  		INTEGER(C_INT) FUNCTION colamd_recommended_c(a, b, c) BIND(C,NAME='colamd_recommended')
            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT
            IMPLICIT NONE
            INTEGER(C_INT), VALUE :: a, b, c
        END FUNCTION colamd_recommended_c
  		INTEGER(C_INT) FUNCTION colamd_c(n_row, n_col, Alen, A, p, knobs, stats) BIND(C,NAME='colamd')
            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
            IMPLICIT NONE
            INTEGER(C_INT), value :: n_row, n_col, Alen
            INTEGER(C_INT) :: A(*), p(*)
            REAL(C_DOUBLE) :: knobs(*)
            INTEGER(C_INT) :: stats(*)
        END FUNCTION colamd_c
  		SUBROUTINE colamd_set_defaults_c(knobs) BIND(C,NAME='colamd_set_defaults')
            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_DOUBLE
            IMPLICIT NONE
            REAL(C_DOUBLE) :: knobs(*)
        END SUBROUTINE colamd_set_defaults_c
  		SUBROUTINE colamd_report_c(stats) BIND(C,NAME='colamd_report')
            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT
            IMPLICIT NONE
            INTEGER(C_INT) :: stats(*)
        END SUBROUTINE colamd_report_c
  		INTEGER(C_INT) FUNCTION symamd_c(n, A, p, perm, knobs, stats) BIND(C,NAME='symamd_wrapper')
            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
            IMPLICIT NONE
            INTEGER(C_INT), value :: n
            INTEGER(C_INT) :: A(*), p(*), perm(*)
            REAL(C_DOUBLE) :: knobs(*)
            INTEGER(C_INT) :: stats(*)
        END FUNCTION symamd_c
  		SUBROUTINE symamd_report_c(stats) BIND(C,NAME='symamd_report')
            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT
            IMPLICIT NONE
            INTEGER(C_INT) :: stats(*)
        END SUBROUTINE symamd_report_c
    END INTERFACE

    ALEN = colamd_recommended_c(A_NNZ,A_NCOL,A_NROW)

    ! ALLOCATE (use C indexing start at zero)
    ALLOCATE(A(0:ALEN-1),p(0:A_NCOL),B(0:3),q(0:5))
    A = 0; p = 0; B = 0; q = 0
    A(0:10) = [0,1,4,2,4,0,1,2,3,1,3]
    p(0:4)  = [0,3,5,9,A_NNZ]
    B(0:3)  = [1,2,3,4]
    q(0:5)  = [0,1,3,3,4,B_NNZ]
    ALLOCATE(perm(0:B_N),stats(COLAMD_STATS),knobs(COLAMD_KNOBS))

    ! DUMP THE INPUT MATRIX A
    WRITE(6,'(A,I3,A,I3,A)') "COLAMD",A_NROW,"-BY-",A_NCOL," INPUT_MATRIX:"
    DO col = 0,A_NCOL-1
        length = p(col+1) - p(col) 
        WRITE(6,'(A,I3,A,I3,A)') "COLUMN",col,", WITH ",length," ENTRIES:"
    	DO pp = p(col),p(col+1)-1
            row = A(pp)
            WRITE(6,'(A,I3)') "    ROW ",row
        END DO
    END DO

    ! Default the knobs
    CALL colamd_set_defaults_c(knobs)

    ! ORDER THE MATRIX
    OK = colamd_c(A_NROW,A_NCOL,ALEN,A,p,knobs,stats)
    CALL colamd_report_c(stats)
    IF (OK .ne. 1) THEN
        WRITE(6,'(A,I3)') "COLAMD_ERROR: ",ok
        STOP 
    END IF

    ! Print the column ordering
    WRITE(6,'(A)') 'COLAMD COLUMN ORDERING:'
    WRITE(6,'(A,I3)') '1ST COLUMN: ',p(0)
    WRITE(6,'(A,I3)') '2ND COLUMN: ',p(1)
    WRITE(6,'(A,I3)') '3RD COLUMN: ',p(2)
    WRITE(6,'(A,I3)') '4TH COLUMN: ',p(3)

    ! DUMP THE STRICKLY LOWER TRIANGULAR PART OF THE SYMMETRIC INPUT MATRIX B
    WRITE(6,*)
    WRITE(6,*)
    WRITE(6,'(A,I3,A,I3,A)') 'SYMAMD ',B_N,'-by-',B_N,' INPUT MATRIX:'
    WRITE(6,'(A)') 'ENTRIES IN STRICTLY LOWER TRIANGULAR PART:'
    DO col = 0,B_N-1
        length = q(col+1) - q(col) 
        WRITE(6,'(A,I3,A,I3,A)') "COLUMN",col,", WITH ",length," ENTRIES:"
    	DO pp = q(col),q(col+1)-1
            row = B(pp)
            WRITE(6,'(A,I3)') "    ROW ",row
        END DO
    END DO

    ! Default the knobs
    CALL colamd_set_defaults_c(knobs)

    ! ORDER THE MATRIX
    OK = symamd_c(B_N,B,q,perm,knobs,stats)
    CALL symamd_report_c(stats)
    IF (OK .ne. 1) THEN
        WRITE(6,'(A,I3)') "SYMAMD_ERROR: ",ok
        STOP 
    END IF

    ! Print the symmetric ordering
    WRITE(6,'(A)') 'SYMAMD COLUMN ORDERING:'
    WRITE(6,'(A,I3)') '1ST COLUMN: ',perm(0)
    WRITE(6,'(A,I3)') '2ND COLUMN: ',perm(1)
    WRITE(6,'(A,I3)') '3RD COLUMN: ',perm(2)
    WRITE(6,'(A,I3)') '4TH COLUMN: ',perm(3)
    WRITE(6,'(A,I3)') '5TH COLUMN: ',perm(4)


    DEALLOCATE(A,p,B,q,perm,stats,knobs)
END PROGRAM ALGO836_TEST